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Abstract. - Within the framework of the mean-field-hydrodynamic model of a degenerate Fermi 
gas (DFG), we study, by means of numerical methods and variational approximation (VA), the 
formation of fundamental gap solitons (FGSs) in a DFG (or in a BCS superfluid generated by 
weak interaction between spin-up and spin-down fermions) , which is trapped in a periodic optical- 
lattice (OL) potential. An effectively one-dimensional (ID) configuration is considered, assuming 
strong transverse confinement; in parallel, a proper ID model of the DFG (which amounts to the 
known quintic equation for the Tonks- Girardeau gas in the OL) is considered too. The FGSs 
found in the first two bandgaps of the OL-induced spectrum ( unless they are very close to edges 
of the gaps) feature a (tightly-bound) shape, being essentially confined to a single cell of the OL. 
In the second bandgap, we also find antisymmetric tightly-bound subfundamental solitons (SFSs), 
with zero at the midpoint. The SFSs are also confined to a single cell of the OL, but, unlike the 
FGSs, they are unstable. The predicted solitons, consisting of ~ 10 4 — 10 5 atoms, can be created 
by available experimental techniques in the DFG of 6 Li atoms. 
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Introduction: Matter-wave solitons were created as lo- 
calized nonlinear excitations in Bosc-Einstcin condensates 
(BECs) of attractively interacting 7 Li [1] and 85 Rb [2] 
atoms loaded in a cigar-shaped trap. In both cases, the 
interaction was switched from repulsion to attraction by 
applying magnetic field near the Feshbach resonance [3,4]). 
In the absence of axial confinement, the solitons can prop- 
agate freely in the axial direction. This was followed by 
the creation of gap solitons, GSs (formed by a few hun- 
dred atoms of 87 Rb) in a self-repulsive BEC loaded in 
a cigar-shaped trapped combined with an optical lattice 
(OL) [5], which was created as the interference pattern 
by counter-propagating laser beams (see also review [6]). 
Theoretical description of the dilute BEC relies upon the 
Gross-Pitaevskii equation (GPE), which provides for a re- 
markably accurate description of various matter- wave pat- 
terns, including solitons [7]. In particular, GSs exist at 
values of the chemical potential that fall in finite gaps of 
the band spectrum of the linear problem, induced by the 
periodic OL potential. In that case, the system possesses 
a negative effective mass, which allows the formation of 
solitons in balance with the repulsive nonlinearity [8]. 



The quest for solitons in degenerate Fermi gases 
(DFGs), which are also available to experiments [9], is a 
new challenge. Quasi-soliton excitations in DFGs formed 
by nonintcracting fermions were analyzed in refs. [10]. 
Solitons were predicted in fermion-boson mixtures [11, 12] 
which feature strong fermion-boson attraction. Currently, 
6 Li- 7 Li [13], 6 Li- 23 Na [14], and 40 K- 87 Rb [15] mixtures of 
boson and fermion atoms are available, as well as binary 
DFGs, such as mixtures of two different spin states in 
40 K [16] and 6 Li [17] gases. Accordingly, solitons in a bi- 
nary fermion gas with attraction between the two species 
have been predicted [18]. 

In this Letter, we aim at predicting solitons of the gap 
type in a DFG trapped in the OL, using a mean- field 
hydrodynamic (MFHD) model of the DFG, which actu- 
ally stems from an approximation of the Thomas-Fermi 
type [19,20]. While the MFHD equations do not grasp 
effects of the multi-particle coherence, they provide a rea- 
sonably accurate description of various macroscopic pat- 
terns in DFGs, including the formation of solitons [21] and 
miscibility-immiscibility transition in binary gases [22]. As 
explained below, essentially the same equations may also 
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apply to a different setting, viz., a BCS (Bardeen-Cooper- 
Schrieffer) superfluid formed in the gas of fermions due to 
a weak attraction between atoms with opposite polariza- 
tions of the spin; the attraction may be induced by means 
of the Fcshbach-rcsonancc technique too. 

We consider three versions of the ID MFHD equation 
for the fermion wave function. The most fundamental 
one is derived from the underlying 3D Fermi distribu- 
tion, assuming a quasi- ID (cigar-shaped) trap with the 
radius of the transverse confinement, a±, larger than the 
de Broglie wavelength, Xf, of atoms on the Fermi sur- 
face. In this case, the MFHD equation for wave function 
-0 contains the self-repulsive nonlinear term |?/>| 4 ' /3 '0. In 
the opposite case of an extremely tight transverse confine- 
ment, with a± <C Xp, the underlying Fermi distribution is 
one-dimensional, leading to an equation with the repulsive 
quintic term, | - i/ , | 4 V'i the same as in the version of the GPE 
derived in ref. [23] for the Tonks-Girardeau gas. Solitons 
in the ID equation combining the quintic term and the 
OL potential were recently considered in refs. [24]. One 
may also consider an anisotropic trap, with aj_ > Xp in 
one transverse direction and a± <C Xp in the other. Then, 
the underlying 2D Fermi distribution gives rise to an equa- 
tion with the self-repulsive cubic term, |V'| 2 '0j formally the 
same as in the ordinary BEC. 

GS solutions to the GPE with the repulsive nonlinearity 
are usually assumed to be loosely bound ones, featuring 
weak localization and wavy tails, see, e.g., refs. [24]. We 
are interested in a possibility to predict GSs of a different 
kind, tightly-bound ones, localized practically in a single 
cell of the OL (cf. refs. [25], [26], and [27], where both 
loosely and tightly bound GSs were investigated in BEC 
models). We demonstrate that this is the case indeed, 
except when the chemical potential is very close to edges 
of the bandgap (in that case, GSs start to develop wavy 
tails, signalizing a transition to delocalized states in the 
Bloch bands). We report families of symmetric (even) 
fundamental gap solitons (FGSs), originating in the first 
bandgap and continuing into the second gap, which feature 
a waveform with a single maximum, and antisymmetric 
(odd) subfundamental gap solitons (SFSs, the name bor- 
rowed from ref. [25]), that emerge in the second bandgap, 
with two maxima in the density profile and a zero between 
them, all squeezed into a single OL cell. 

Below, we study these solitons by means of numerical 
simulations and variational approximation (VA) [28] . The 
VA, based on the Gaussian ansatz, predicts the FGS fami- 
lies, in terms of the dependence between the effective non- 
linearity and chemical potential, in the first two bandgaps 
with a surprisingly high accuracy, if compared to numer- 
ical results. SFSs arc predicted too, although less accu- 
rately, in the second bandgap, by means of another (an- 
tisymmetric) ansatz. We also report examples of stable 
symmetric and antisymmetric bound states of the FGSs. 

The models We derive the MFHD equation for the 
fermion wave function, as an extension of the static 
Thomas-Fcrmi distribution of the atomic density, n = 



|^| 2 , which is [19,20] 

V{r) + Vg = M, (1) 

with fj, the chemical potential, V(r) the external trap- 
ping potential, and V^g an effective nonlinear potential 
in the P-dimcnsional space accounting for the Fermi pres- 
sure. For T> = 3D, 2D, and ID, the role of \Q S is ac- 
tually played by the local Fermi energy, e-p, expressed 
in terms of local atomic density up, which is defined 
as per the dimension [20]. This yields (2m/ ti 2 ) = 
(67r 2 n3D) 2 / 3 , 4-7m2D, and (imm) 2 , respectively, with m the 
atom mass. The MFHD equation for stationary problems 
is derived by treating relation (1) as one which defines 
the effective Thomas-Fermi Hamiltonian in terms of V^. 
Augmenting the Hamiltonian with the kinetic-energy term 
and applying it to the wave function, one arrives at the 
stationary equation [20], 

A dynamic version of eq. (2) is [20] 

h 2 Bth 

~2^^ J + [v{r) + v ^ = lh lTv (3) 

wither, t) = e _vt /' i ^ r (r), although this generalization is 
formal, as the above-mention derivation docs not define 
the phase of the fermion wave function. Nevertheless, the 
dynamical equation, if it is treated as a phenomenologi- 
cally postulated one, may yield reasonable predictions for 
stability of static MFHD states in the DFG [12,18,20]. 

As mentioned above, dynamical equation (3) can also 
be derived in a different physical context, for a BCS su- 
perfluid formed by Cooper pairs of fermions with opposite 
polarizations of the spin. Indeed, using the known ex- 
pression for the local energy density of the superfluid in 
3D [29], £30 = (3/5)n 3 D£F (recall e-p is the local Fermi 
energy), and deriving the equation for the wave function 
of the Cooper pairs (alias the complex order parameter, 
in terms of the Ginzburg-Landau approach) from the La- 
grangian density which includes term £sd, one will end 
up with a 3D equation that differs from the 3D version of 
eq. (3) only by a numerical factor in front of V^g. In a 
similar fashion, using the local density of the BCS super- 
fluid in ID [30] and 2D settings, £id = (l/3)n 3 D£F and 
<?2D = (1/2)ti2d£f, one can derive the respective ID and 
2D equations, which differ from their counterparts (7) and 
(6) written below only by numerical coefficients in front 
of the nonlinear terms. 

In the case of the ordinary cigar-shaped waveg- 
uide with strong harmonic confinement in the trans- 
verse plane, eq. (3) and its static version, eq. 
(2), can be reduced to the ID form by means 
of the well-known substitution [31], tp(x,y, z,t) = 
</>(x, t) exp [— iu±t — (y 2 + z 2 ) / (2aj_)] ,with the second 
multiplier representing the ground state of the transverse 
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harmonic oscillator, uj± is the respective frequency, and 
aj_ = \J K/Jmcj±). The substitution of this ansatz in eq. 
(3) for T> = 3D and averaging in the transverse plane yield 
the effective ID equation, 



ifuj) t = - — 



2m 



3ft 2 
10m 



6tt 2 \<f>\ 



2/3 



'47T 
[ — X 



(4) 

where the OL potential with strength e and period A/2 
is introduced. Equation (4) is further rescaled by defin- 
ing 4> = ^/2N/XaJ 1 4> 1 t = \ 2 mj (47r 2 ft)f, x = A.t/(2tt), 
Vq = em)? J (2irK) 2 , where N is the number of atoms. The 
result is 



i<t>t 



(l/2)<t> xx + ff3 D |0| 4/3 <t>-V cos (2a:) (5) 



(tildes are dropped here), with the wave function sub- 
ject to normalization Af = \(f>(x)\ 2 dx = 1, and ef- 
fective strength of the nonlincarity and potential defined 

as 53D = (3/10) [3A^A 2 /(27ra 2 L )] 2/3 (subscript 3D refers to 
the derivation of the equation from the 3D Fermi distri- 
bution), Vo = me (\/2irh) 2 . Similarly, the following nor- 
malized ID equations can be derived from the underlying 
2D and ID Fermi distributions: 



i<Pt = 



g-2D \<P\ 2 <P- Vo COS (2x) </>, (6) 
■Sid \4>t<t>- Vb cos (2z) </>, (7) 



with 52D = XN/ (y/ra±) , gm = tt 2 N 2 /2. 

Stationary solutions to eqs. (5), (6), and (7) are looked 
for in the usual form, <j>(x,t) = e~~ lfJ - t u(x), with real func- 
tion u obeying the equation 



jiu + (l/2)u" - gu** + V cos (2x) u = 0, 



(8) 



where = 7/3, 3, 5, respectively, for V = 3D, 2D, ID. 
In fact, the static ID equations (8) for the real wave func- 
tion have the straightforward physical meaning, within the 
framework of the MFHD description, while their dynamic 
counterparts for the complex wave function are more for- 
mal ones, as mentioned above. Nevertheless, the dynam- 
ical equations are quite useful, as their direct simulations 
converge to real stationary states that represent numeri- 
cally exact solutions to the static equations (see below). 
Besides that, the dynamical equations make sense (simi- 
lar to the time-dependent Ginzburg-Landau equations) if 
applied, as mentioned above, to the BCS superfiuid. Our 
main objective is to find a family of FGS solutions to eq. 
(8) with /j, falling in the first two finite bandgaps of the 
linear spectrum. 

Variational approximation Equation (8) and the above 
normalization condition, Af = 1, can be derived as the 
variational equations [28], SL/6u = dL/d/j, = 0, from the 
Lagrangian, 



2 1 i /x2 2ff 
2 K-D + 



+Vb cos(2a;) ■ u 2 ] dx — fx. 



(9) 



To apply the VA, we use a simple Gaussian ansatz [28] , 



u{x) = Tr-^V-WyVKcxp [_a;2y ( 2 VF 2 )] 



(10) 



where variational parameters are the soliton's norm and 
width, Af and W (in addition to fi). The substitution of 
this ansatz in Lagrangian (9) yields 



L = n {M - 1) - 



Af 



V AfNe 



-w* 



AW 2 

2 3 / 2 g W(^+i)/2 

n Kv/2 (H P + 1) 3/2 W^-l)/2- 



(11) 



The corresponding variational equations, dLjd[i = 
dL/dW = dL/dAf = 0, lead to Af = 1 (as expected), 
and 



1 + 2^(H P -l) g _, p)/2 = 4g _ 

7r(K»-l)/4 (Np + 1) 3/2 



ir- 



1 



A* = 



V2g 



AW 2 



-V n e- 



(12) 
(13) 



Then, eqs. (12) and (13) were solved numerically. 

Results Simulations of eqs. (5), (6), and (7) were car- 
ried out by dint of a real-time integration method based 
on the Crank-Nicholson discretization scheme, as elabo- 
rated in rcf. [32]. The equations were discretized using 
time and space steps 0.0005 and 0.025, respectively, in 
domain —20 < x < 20. To find FGSs, the simulations 
started with an initial configuration chosen as the ground 
state, <f>(x) — (V2c/ 7r) 1 / 4 exp[— x 2 y/c/2], of the linear os- 
cillator with potential cx 2 (with c 3> 1, typically 5 to 100). 
The OL potential was slowly introduced in the course of 
the simulations. The strong harmonic potential squeezes 
the localized state into a single cell of the OL. After ob- 
taining a stationary bound state in the nonlinear equation 
containing the combined OL and harmonic potential, the 
latter one was slowly switched off. The simulations were 
run until a well-defined stationary shape of the soliton was 
established. A similar approach to the creation of GSs in 
BEC was proposed in ref. [33]. To generate SFS solutions 
(see below), the initial state was taken as the first excited 
state of the above-mentioned harmonic potential, instead 
of its ground state. We present results for a characteristic 
value of the OL strength, Vb = 5. 

Families of FGS solutions are characterized by the corre- 
sponding dependences, <7x>(/x). In fig- 1, they are displayed 
as found from the numerical solution of the ID, 2D, and 
3D models, along with the same dependences as predicted, 
for the three models, by the VA. 

In fig. 2 we display typical profiles of the solitons. It is 
obvious that the FGSs, unless taken too close to bandgap 
edges, arc compact objects, primarily trapped in a single 
cell of the OL. There also exist more loosely organized 
localized states, which occupy several cells [24], that we 
do not consider here. Figure 2 demonstrates that the VA 
produces a very good fit not only to the g(fjf) plots for the 
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Fig. 1: Numerical (continuous curve) and variational (chain 
of symbols) results for nonlinearity g versus chemical poten- 
tial fj,, for fundamental gap solitons in the first and second 
finite bandgaps of periodic potential — Vb cos(2a;), (all results 
reported in this Letter with V"o = 5,) in the models based on 
eqs. (7), (6), and (5) ("ID", "2D", and "3D", respectively, 
which refers to the dimension of the underlying Fermi distri- 
bution, from which the equation is derived) . The shaded areas 
represent the Bloch bands. 



entire soliton families (see fig. 1), but also to FGS profiles 
in the two lowest finite gaps, except very close to edges of 
the gaps, where the solitons develop undulate tails, that 
simple ansatz (10) is unable to reproduce [see panels in 
fig. 2 pertaining to g = 0.01]. In principle, GSs with 
conspicuous tails may be approximated by a more com- 
plex ansatz, which combines the Gaussian and function 
cos a;; however, variational equations generated by such 
an ansatz are very cumbersome, and, in the end, the ex- 
tended ansatz does not approximate the entire tail, but 
only its secondary peaks which are closest to the central 
one [27]. For this reason, we do not try to apply that 
ansatz here. 

It is relevant to mention that as the chemical potential 
moves closer to the Bloch band the wave form of the GSs 
develop a structure in space similar to Wannier functions, 
which are spatially localized linear combinations of Bloch 
functions. Figure 2 (a) shows the wave form near the 
edges of the gap. The undulate tails in this figure signal 
the proximity to a Bloch band, where the wave form will 
turn into a periodic structure. However, as said above, 
such wave forms cannot be fitted to Gaussian shapes and 
we do not intend to study them in detail here. 

Within the framework of eqs. (5) - (7), the stability 
of the FGSs was tested by subjecting them to relatively 
strong initial perturbations (as said above, the dynamical 
equations have direct meaning in the application to the 
BCS supcrfluid). It has been found that the entire fam- 
ilies of these solutions are stable, in all the three above- 
mentioned models, see a typical example in Fig. 3. In 
the case displayed in this figure, after the stationary FGS 
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Fig. 2: Typical shapes of the fundamental gap solitons in the 
model derived from the ID Fermi distribution, Eq. (7), for 
(a) g m = 0.01, (/i = -2.892), (b) g m = 20 (/x = 2.427), 
and from the 3D distribution, Eq. (5), for (c) gsu = 0.01, 
(At = -2.891), and (d) g 3 u =9,(n = 3.284 ). Shapes produced 
by the numerical solution and variational approximation are 
labeled as "num" and "var". The panels pertaining to gsu — 
gin =0.01 display the solitons very close to the left edge of the 
first bandgap, while the panels for gio = 20 and gw = 9 give 
typical examples found inside the second bandgap (inside the 
first gap, the shape of the fundamental solitons is quite similar) . 
The thin sinusoidal line (in this and following figures) depicts 
the profile of the underlying OL potential. 



of the model based on the 3D Fermi distribution was ob- 
tained (for #3D = 9), we replaced, at t — 20, the stationary 
wave form <f>(x) by a perturbed one, 1.2<f>(x), and contin- 
ued the simulation. The resultant solution demonstrates 
persistent oscillations, and no sign of degradation. 

The usual GPE, with the repulsive cubic nonlinearity 
and OL potential, i.e., eq. (6), gives rise to antisymmet- 
ric SFSs, the family of which starts at the left edge of 
the second finite bandgap [25]. A characteristic feature 
of the SFSs is that two maxima of the density, \4>(x)\ , 
are located inside a single site of the periodic potential 
(i.e., this antisymmetric soliton as a whole is confined to 
a single site). To look for SFSs in the present models, we 
applied the VA based on the modified Gaussian ansatz (cf. 
cq. (10 )), 



u(x) 



xexp {—x 2 / (2W 2 )) , (14) 



where N and W have the same meaning as in cq. (10). 
The substitution of ansatz (14) in Lagrangian (9) yields 



L = /* CAT — 1) 



3A/~ 
4W 2 



+ V Ne~ w {1-2W 2 ) 



(15) 



7r(»o+i)/4 (Hp + i)^+ 4 )/ 2 M/(^~i)/2 ; 
and the variational equations following from here are M = 



7-\_/l 
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Fig. 3: An example of stable evolution of a fundamental gap 
soliton solution in the 3D model [eq. (5)], for g3u — 9 of Fig. 2 
(d). To introduce a perturbation, at t = 20, the wave form was 
suddenly multiplying by 1.2: cj>(x,t) — > 1.2 x cj>(x,t). The so 
perturbed soliton (as well as ones subjected to different arbi- 
trary perturbations) remained stable as long as the simulations 
were run. 




1 and 



Hp - 1 2 {m ' D+5)/2 gW {5 ~^ T,)/2 
3 ^+D/4 (Hp + 1)< 4 +^/ 2 



= (4/3)V W 



(3-2W 2 ), 



r(Hp/2 + i) 

(16) 



2(2K D +3)/2 



m v /2+i) 



AW 2 



tt^o+i)/* (Hp + i)( 2 +^)/ 2 IF^- 1 )/ 2 



(17) 



Numerical solution of eqs. (16) and (17) gives rise to 
SFS families in the second bandgap. The comparison of 
the corresponding characteristics g(p-) with their numeri- 
cally found counterparts is presented in fig. 4(a). The ac- 
curacy of the VA predictions for the SFS families is worse 
than for the FGS solutions; nevertheless, the prediction is 
qualitatively correct. 

Although the SFS families were found from direct sim- 
ulations of eqs. (5), (6), and (7), hence they are, at least, 
quasi-stable solutions, longer simulations reveal their in- 
stability (not shown here). Similar to what was found in 
ref. [25] in the GPE with the cubic nonlinearity, which 
coincides with eq. (6), in all the three models considered 
here, the instability tends to transform the SFS (double- 
humped) solitons into a stable FGS (single- humped). 

In addition to the SFSs, symmetric and antisymmetric 
( "twisted" , in the latter case) bound states of FGSs were 
found too. Typical examples of such stable bound states in 
the 3D and ID models are shown in figs. 5(a) and (b), re- 
spectively. The antisymmetric state is formally similar to 
the SFS, but a difference is that the density maxima of the 
bound state are located in different OL sites. Direct sim- 
ulations (not shown here) clearly demonstrate that both 



Fig. 4: (a) The same as in fig. 1, but for families of subfun- 
damental gap solitons. (b) A typical example of the subfunda- 
mental soliton, obtained from the numerical solution of eq. (5) 
(for (/3d = 3), and its counterpart predicted by the VA. 



symmetric and antisymmetric bound states are stable one 
(on the contrary to the SFS solutions). 

Conclusion In the three MFHD models, based on eqs. 
(7), (6), and (5), we have considered, by means of nu- 
merical simulations and VA (variational approximation), 
the generation of stable fundamental and unstable sub- 
fundamental gap solitons (FGSs and SFSs, respectively) 
in the quasi-lD degenerate Fermi gas loaded in a periodic 
OL potential. The same models apply to the BCS super- 
fluid trapped in the OL. In most cases, both the FGSs 
and SFSs are confined, essentially, to a single cell of the 
lattice. The VA provides a very good fit to the FGSs, 
and a qualitatively correct approximation for the SFSs. 
Stable symmetric and antisymmetric bound states of the 
FGSs were also found. It was demonstrated that trans- 
port of the FGSs, without much distortion, is possible by 
a slowly moving OL. 

Experimental realization of the FGSs in degenerate 
Fermi gases seems quite feasible. To this end, the gas 
should be loaded in a cigar-shaped trap combined with 
the periodic axial potential, as was done to create gap 
solitons in the BEC [5,6]. The experiment may be started 
with an additional strong parabolic potential acting in the 
axial direction, to prepare a strongly localized state, that 
may have a good chance to self-trap into a tightly bound 
soliton while the extra potential is gradually switched off. 
An estimate for the 6 Li gas, with typical values of phys- 
ical parameters, predicts that the FGSs based on the 3D 
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Fig. 5: Examples of stable symmetric (a) and antisymmetric 
(b) bound states of fundamental gap solitons in the models 
based on the 3D [Eq. (5)] and the ID [Eq. (7)] Fermi distri- 
butions, respectively, for (a) g^ — 10 and (b) gio = 6. 



Fermi distribution, i.e., obeying eq. (5), may be created 
with 10 4 — 10 5 atoms per soliton. 
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